

package scalaSci

import scalaSci.math.LinearAlgebra.LinearAlgebra
import scalaSci.math.array.DoubleArray
import java.util._
import Jama._


// this class provides Matrix operations, using zero-indexed values
class Mat( n: Int, m: Int) {    // the default constructor allocates a double array of size n x m
   private var  N: Int=n
   private var  M: Int=m
   private  var  v = new Array[Array[Double]](N, M)  // i.e. v(1:n, 1:m)

  def  setv(values: Array[Array[Double]], n: Int, m: Int) = { v = values; N = n; M = m }
   
  def size = {  (N, M) }  // the size of the Matrix

  def length =  { N }

  // getters for size
  def Nrows = { N }
  def Ncols = { M }

  def getv = {
    v
}
    // construct a Matrix from a zero-indexed double [][] array
def this(vals: Array[Array[Double]]) = {
     this(vals.length, vals(0).length)  // allocate memory with the default constructor
    var i=0; var j=0
    while (i<vals.length) {
        j=0
        while (j<vals(0).length) {
       v(i)(j) = vals(i)(j)
       j+=1
        }
        i+=1
     }
 }



def this(N: Int, M: Int, df: Double) = {   // initializes to a default value
     this(N, M)  // allocate memory with the default constructor
    var i=0; var j=0
    while (i<v.length) {
        j=0
        while (j<v(1).length) {
       v(i)(j) = df
       j+=1
        }
        i+=1
     }
 }

    // construct a Matrix from an one-indexed double [][] array
def this(vals: Array[Array[Double]], flag: Boolean) = {
    this(1,1)
    v = vals
    N = vals.length
    M = vals(0).length
 }

  // construct a Matrix from a one-indexed double [] array
def this( vals: Array[Double]) = {
    this(1, vals.length)   // keep the array as a row of the Matrix
    var i=0
    while (i<vals.length)  {
       v(0)(i) = vals(i)
       i+=1
   }
  }

    // clone a Matrix from another Matrix
def this(a: Mat) = {
     this(a.N, a.M)  // allocate memory with the default constructor
    var i=0; var j=0
    while (i<a.N) {
        j=0
        while (j<a.M) {
       v(i)(j) = a.v(i)(j)
       j+=1
        }
        i+=1
     }
 }

  // apply the function f to all the elements of the Mat and return the results with a new Mat
 def  map( f: (Double => Double)): Mat = {
   var mres = new Mat(N, M)
   
    for (r<-0 until N)
      for (c<-0 until M) 
       mres(r, c) = f(this(r, c) )
   
   mres
 }

def  print: String = {
  var digitFormat = scalaExec.Interpreter.GlobalValues.fmtMatrix
    var nelems = N
    var melems = M

    var s:String =" \n"
     var i=0; var j=0;
    while (i<nelems) {
        j=0
        while (j<melems) {
       s = s+digitFormat.format(v(i)(j))+"  "
       j+=1
        }
        s += "\n";
     i+=1
    }
  s
 }

override  def  toString(): String = {
  var digitFormat = scalaExec.Interpreter.GlobalValues.fmtMatrix
    var nelems = N
    var melems = M
    var truncated = false

    var truncationIndicator = ""
    if  (Mat.mxElemsToDisplay < N) {
        nelems = Mat.mxElemsToDisplay
        truncationIndicator = " ... "
        truncated = true
      }
     if  (Mat.mxElemsToDisplay < M) {
        melems = Mat.mxElemsToDisplay
        truncationIndicator = " ... "
      }
     var s:String =" \n"
     var i=0; var j=0;
    while (i<nelems) {
        j=0
        while (j<melems) {
       s = s+digitFormat.format(v(i)(j))+"  "
       j+=1
        }
      s += (truncationIndicator+"\n")
     i+=1
    }
   if (truncated)     // handle last line
    s+= "........................................................................................................................................."

     return s
 }

def p = print   // short "print" method 

    // ROW - APPEND MAT TO MAT
// append rowwise the Matrix M
def  ::(rowsToAppend: Mat): Mat = {
    if (rowsToAppend.M != M )   // incompatible number of columns
      return this
    // create a new extended matrix to have also the added rows
    var  exrows = N+rowsToAppend.N   // new number of rows
    var  res = new Mat(exrows, M)
    // copy "this" Matrix
    for (r<-0 to N-1)
       for (c<-0 to M-1)
         res.v(r)(c) = v(r)(c)
    for (r<-0 to rowsToAppend.N-1)
       for (c<-0 to rowsToAppend.M-1)
         res.v(N+r)(c) = rowsToAppend.v(r)(c)
    res
}


    // ROW - PREPEND MAT TO MAT
// prepend rowwise the Matrix M
def  :::(rowsToAppend: Mat): Mat = {
    if (rowsToAppend.M != M )   // incompatible number of columns
      return this
    // create a new extended matrix to have also the added rows
    var  exrows = N+rowsToAppend.N   // new number of rows
    var  res = new Mat(exrows, M)
    // copy "this" Matrix
    for (r<-0 to rowsToAppend.N-1)
       for (c<-0 to rowsToAppend.M-1)
         res.v(r)(c) = rowsToAppend.v(r)(c)
    for (r<-0 to N-1)
       for (c<-0 to M-1)
         res.v(rowsToAppend.N+r)(c) = v(r)(c)

    res
}


    // COL - APPEND MAT TO MAT
// append columnwise the Matrix M
def  >::(colsToAppend: Mat): Mat = {
    if (colsToAppend.N != N )   // incompatible number of rows
      return this
    // create a new extended matrix to have also the added columns
    var  excols = M+colsToAppend.M   // new number of columns
    var  res = new Mat(N, excols)
    // copy "this" Matrix
    for (r<-0 to N-1)
       for (c<-0 to M-1)
         res.v(r)(c) = v(r)(c)
    for (r<-0 to colsToAppend.N-1)
       for (c<-0 to colsToAppend.M-1)
         res.v(r)(M+c) = colsToAppend.v(r)(c)
    res
}


    // COL - PREPEND MAT TO MAT
// prepend columnwise the Matrix M
def  >:::(colsToAppend: Mat): Mat = {
    if (colsToAppend.N != N )   // incompatible number of rows
      return this
    // create a new extended matrix to have also the added columns
    var  excols = M+colsToAppend.M   // new number of columns
    var  res = new Mat(N, excols)
    // copy "this" Matrix
    for (r<-0 to colsToAppend.N-1)
       for (c<-0 to colsToAppend.M-1)
         res.v(r)(c) = colsToAppend.v(r)(c)
    for (r<-0 to N-1)
       for (c<-0 to M-1)
         res.v(r)(colsToAppend.M+c) = v(r)(c)

    res
}


// indexes the corresponding Matrix element without updating automatically sizes for fast access
  def apply(n: Int, m: Int): Double = {
        v(n)(m)

 }


  // indexes the corresponding Matrix element, resizing automatically 
  def apply(n: Int, m: Int, resize: Boolean): Double = {
      if (n>=N || m >= M)  {  // update size
          var newN = (Mat.matResizeFactor*n).asInstanceOf[Int]
          var newM = (Mat.matResizeFactor*m).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // update size

        v(n)(m)

 }



// extracts a submatrix specifying rows only, take all columns, e.g. m(2, 3, ':') corresponds to Matlab's m(2:3, :)'
// m(low:high,:) is implemented with m(low, high, dummyChar). if low>high then rows are returned in reverse
 def apply(rowL: Int, rowH: Int, allColsChar: Char): Mat = {
   var rowStart = rowL; var rowEnd=rowH;
   var colStart = 0;     var colEnd =  M-1;   // all columns
   var colNum = M
   var colInc = 1

if (rowStart <= rowEnd) {   // positive increment
    var rowInc = 1
    if (rowEnd == -1) { rowEnd = N-1 }  // if -1 is specified take all the rows
    var rowNum = rowEnd-rowStart+1
    var subMatr = new Mat(rowNum, colNum)   // create a Mat to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart
    var rowIdx =0; var colIdx = 0  // indexes at the new Matrix
    while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd )   { 
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       } // crow <= rowEnd
subMatr   // return the submatrix

} // rowStart <= rowEnd
else { // rowStart > rowEnd
    var rowInc = -1
    var rowNum = rowStart-rowEnd+1
    var subMatr = new Mat(rowNum, colNum)   // create a Mat to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row at the source matrix
    var ccol = colStart
    var rowIdx =0; var colIdx = 0  // indexes at the new Mat
    while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       }

subMatr   // return the submatrix

} // rowStart > rowEnd

}

// extracts a submatrix specifying rows only, take all columns, e.g. m(2, 4, 12, ':') corresponds to Matlab's m(2:4:12, :)'
def apply(rowL: Int, rowInc: Int, rowH: Int, allColsChar: Char): Mat = {
    var rowStart = rowL;     var rowEnd =  rowH;
    var colStart = 0;  var colEnd = M-1;   // all columns
    var colNum = M
    var colInc = 1

  if (rowInc > 0) { // positive increment
    if (rowEnd == -1) { rowEnd = N-1 }  // if -1 is specified take all the rows
    var rowNum = Math.floor( (rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
    var colStart =0;     var colEnd =  M-1   // all columns
    var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart
    var rowIdx = 0; var colIdx = 0  // indexes at the new Mat
    while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       }
     subMatr   // return the submatrix
     }  // positive increment
  else  {  //  negative increment
     var rowNum = Math.floor( (rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
     var subMatr = new Mat(rowNum, colNum)  // create a Matrix to keep the extracted range
        // fill the created matrix with values
     var  crow = rowStart   // indexes current row at the source matrix
     var  ccol = colStart
     var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix
     while (crow >= rowEnd)  {
         ccol = colStart;  colIdx = 0
         while (ccol <= colEnd)  {
             subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
             colIdx += 1
             ccol += colInc
         }
         rowIdx += 1
         crow += rowInc
      }
       subMatr  // return the submatrix
     }  // negative increment
}



// extracts a submatrix, e.g. m( ':', 2,  12 ) corresponds to Matlab's m(:, 2:12)'
  def apply(allRowsChar: Char, colLow: Int,  colHigh: Int): Mat = {
   var rowStart = 0;     var rowEnd =  N-1   // all rows
    var colStart = colLow;  var colEnd = colHigh
    var rowInc = 1
    var colInc = 1
    var rowNum = N    // take all the rows

    if  (colStart <= colEnd)   {    // positive increment
        if (colEnd == -1)  { colEnd = M-1 } // if -1 is specified take all the columns
        var colNum = colEnd-colStart+1
        var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += rowInc
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else {  // negative increment
    var colNum = colEnd-colStart+1
    var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
   }



// extracts a submatrix, e.g. m( ':', 2, 3, 12 ) corresponds to Matlab's m(:, 2:3:12)'
  def apply(allRowsChar: Char, colLow: Int, colInc: Int, colHigh: Int): Mat = {
   var rowStart = 0;     var rowEnd =  N-1   // all rows
    var colStart = colLow;  var colEnd = colHigh
    var rowInc=1
    var rowNum = N    // take all the rows

    if  (colStart <= colEnd)   {    // positive increment
        if (colEnd == -1)  { colEnd = M-1 } // if -1 is specified take all the columns
        var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
        var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else {  // negative increment
    var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
    var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
   }


// extracts a submatrix, e.g. m( 2, 3, 12, 4, 2,  8 ) corresponds to Matlab's m(2:3:12, 4:2:8)'
  def apply(rowLow: Int, rowInc: Int, rowHigh: Int, colLow: Int, colInc: Int, colHigh: Int): Mat = {
    var rowStart = rowLow;     var rowEnd =  rowHigh
    var colStart = colLow;  var colEnd = colHigh

        var rowNum = Math.floor((rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
        var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
        var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range

    if  (rowStart <= rowEnd && colStart <= colEnd)   {    // positive increment at rows and columns
        var crow = rowStart  // indexes current row
        var ccol = colStart  // indexes current column
        var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix
            while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else if  (rowStart >= rowEnd && colStart <= colEnd)   {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else if  (rowStart <= rowEnd && colStart >= colEnd)   {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow > rowEnd
 subMatr   // return the submatrix
   }

   }


// extracts a specific row, take all columns, e.g. m(2) corresponds to Matlab's m(2, :)'
  def apply(row: Int): Mat = {
    var colStart = 0;     var colEnd =  M-1;   // all columns
    var rowNum = 1;  var colNum = colEnd-colStart+1;
    var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var ccol = colStart
    while  (ccol <= colEnd)   {
          subMatr.v(0)(ccol) = v(row)(ccol)
          ccol += 1
         }

     subMatr
}

// extracts a specific row, take all columns, e.g. m(2, ':') corresponds to Matlab's m(2, :)'
  def apply(row: Int, allColsChar: Char): Mat = {
      if (allColsChar==':')
          apply(row)
      this
      }

// extracts a specific column, take all rows, e.g. m(':', 2) corresponds to Matlab's m(:,2:)'
  def apply(colonChar:Char, col: Int): Mat = {
    var rowStart = 0;     var rowEnd =  N-1;   // all rows
    var colNum = 1;  var rowNum = rowEnd-rowStart+1;
    var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart
    while  (crow <= rowEnd)   {
          subMatr.v(crow-rowStart)(0) = v(crow)(col)
          crow += 1
         }

     subMatr
}


// extracts a submatrix, e.g. m( 2,  12, 4,   8 ) corresponds to Matlab's m(2:12, 4:8)'
  def apply(rowLow: Int,  rowHigh: Int, colLow: Int, colHigh: Int): Mat = {
    var rowStart = rowLow;     var rowEnd =  rowHigh
    var colStart = colLow;  var colEnd = colHigh
    var rowInc = if (rowHigh > rowLow) 1 else -1
    var colInc = if (colHigh > colLow) 1 else -1

        var rowNum = Math.floor((rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
        var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
        var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range

    if  (rowStart <= rowEnd && colStart <= colEnd)   {    // positive increment at rows and columns
        var crow = rowStart  // indexes current row
        var ccol = colStart  // indexes current column
        var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix
            while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else if  (rowStart >= rowEnd && colStart <= colEnd)   {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else if  (rowStart <= rowEnd && colStart >= colEnd)   {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 0; var colIdx = 0  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 0
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow > rowEnd
 subMatr   // return the submatrix
   }

   }




              // extracts a submatrix, e.g. m(3:2:7, :)
  def apply(rowLow: Int, rowInc: Int, rowHigh: Int): Mat = {
    var rowStart = rowLow;     var rowEnd =  rowHigh;    if (rowEnd < rowStart) { rowStart = rowHigh; rowEnd = rowLow; }
    var colStart = 1;     var colEnd =  M-1;
    var colInc = 1
    var rowNum = Math.floor( (rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
    var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
    var subMatr = new Mat(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var inc = 1
    var ccol = colStart
    var rowIdx = 0; var colIdx = 0;  // indexes at the new Matrix
    while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0;
          while  (ccol <= colEnd)    {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       }
     subMatr
}


// returns the corresponding row of the Mat class as an Array[Double]
def getRow(row: Int): Array[Double] = {
    var colStart = 0;     var colEnd =  M-1;   // all columns
    var colNum = colEnd-colStart+1;
    var rowArray = new Array[Double](colNum)
    for (ccol<-0 to colEnd)
       rowArray(ccol) = v(row)(ccol)
    rowArray
  }



    // updating a single element of the Matrix, without resizing
def  update(n: Int, m: Int, value: Double): Unit = {
   if (n>=N || m >= M)  {  // update size
          var newN = (Mat.matResizeFactor*n).asInstanceOf[Int]
          var newM = (Mat.matResizeFactor*m).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // update size

         v(n)(m) = value;
   }

  def  update(n: Int, m: Int, value: Double, resize: Boolean): Unit = {
   if (n>=N || m >= M)  {  // update size
          var newN = (Mat.matResizeFactor*n).asInstanceOf[Int]
          var newM = (Mat.matResizeFactor*m).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // update size

         v(n)(m) = value;
   }

// update a Matrix subrange by assigning a Matrix, e.g. var m=rand(20, 30);  m( 2, 1, 3, 3, 1, 4) = ones(2,2)
 def update(rlowp:Int, rincp: Int, rhighp: Int, clowp:Int, cincp: Int,  chighp: Int, mr: Mat ): Unit = {
     var rlow=rlowp; var rhigh=rhighp; var rinc = rincp;
     var clow=clowp; var chigh=chighp; var cinc = cincp;
     if (rhigh < rlow)  {
            if (rinc > 0)  {
                println("negative row subrange increment is required")
                this
            }
            var tmp=rhighp; var rhigh=rlowp; rlow = tmp;
            rinc  = -rinc
        }
    if (chigh < clow)  {
            if (cinc > 0)  {
                println("negative column subrange increment is required")
                this
            }
            var tmp=chighp; var chigh=clowp; clow = tmp;
            cinc  = -cinc
        }
      var mrM = mr.M   // length of right-hand side matrix
      var mrN = mr.N
     var rangeLenRow = ((rhigh-rlow)/rinc).asInstanceOf[Int]+1    // row length of target range
     var rangeLenCol = ((chigh-clow)/cinc).asInstanceOf[Int]+1    // col length of target range

     if (mrN != rangeLenRow)  {  // improper sizes of the involved row subranges for assignment
         println("target vector row subrange and right-hand size vector lengths do not agree, i.e. mrN =  "+mrN+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

     if (mrM != rangeLenCol)  {  // improper sizes of the involved column subranges for assignment
         println("target vector column subrange and right-hand size vector lengths do not agree, i.e. mrM =  "+mrM+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

      if (rhigh >= N || chigh >= M)  {   // dynamically increase the size of the matrix when subassigning out of its range
          var newN = (Mat.matResizeFactor*rhigh).asInstanceOf[Int]
          var newM = (Mat.matResizeFactor*chigh).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // dynamically increase the size of the matrix


     // copy the values of the mr
        var rrow=0; var rcol=0; var lrowidx=rlow; var lcolidx = clow
        while (rrow < mr.N) {   // for all rows of the right-hand side matrix
            rcol=0
            lcolidx = clow   // starting column within the "subassigned" matrix
            while (rcol < mr.M)  {   // for all cols of the right-hand side matrix
                v(lrowidx)(lcolidx) = mr(rrow, rcol)
                lcolidx += cinc
                rcol += 1
            }
            lrowidx += rinc
            rrow += 1
          }

        }



// update a Matrix subrange by assigning a Matrix, e.g. var mm = rand(20, 30);  mm(2,3, 3,4) = ones(2,2);
 def update(rlowp:Int, rhighp: Int, clowp:Int, chighp: Int, mr: Mat ): Unit = {
     var rlow=rlowp; var rhigh=rhighp; var rinc = 1;
     var clow=clowp; var chigh=chighp; var cinc = 1;
     if (rhigh < rlow)  {
            if (rinc > 0)  {
                println("negative row subrange increment is required")
                this
            }
            var tmp=rhighp; var rhigh=rlowp; rlow = tmp;
            rinc  = -rinc
        }
    if (chigh < clow)  {
            if (cinc > 0)  {
                println("negative column subrange increment is required")
                this
            }
            var tmp=chighp; var chigh=clowp; clow = tmp;
            cinc  = -cinc
        }
      var mrM = mr.M   // length of right-hand side matrix
      var mrN = mr.N
     var rangeLenRow = ((rhigh-rlow)/rinc).asInstanceOf[Int]+1    // row length of target range
     var rangeLenCol = ((chigh-clow)/cinc).asInstanceOf[Int]+1    // col length of target range

     if (mrN != rangeLenRow)  {  // improper sizes of the involved row subranges for assignment
         println("target vector row subrange and right-hand size vector lengths do not agree, i.e. mrN =  "+mrN+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

     if (mrM != rangeLenCol)  {  // improper sizes of the involved column subranges for assignment
         println("target vector column subrange and right-hand size vector lengths do not agree, i.e. mrM =  "+mrM+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

      if (rhigh >= N || chigh >= M)  {   // dynamically increase the size of the matrix when subassigning out of its range
          var newN = (Mat.matResizeFactor*rhigh).asInstanceOf[Int]
          var newM = (Mat.matResizeFactor*chigh).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // dynamically increase the size of the matrix


     // copy the values of the mr
        var rrow=0; var rcol=0; var lrowidx=rlow; var lcolidx = clow
        while (rrow < mr.N) {   // for all rows of the right-hand side matrix
            rcol=0
            lcolidx = clow   // starting column within the "subassigned" matrix
            while (rcol < mr.M)  {   // for all cols of the right-hand side matrix
                v(lrowidx)(lcolidx) = mr(rrow, rcol)
                lcolidx += cinc
                rcol += 1
            }
            lrowidx += rinc
            rrow += 1
          }

        }




// dot product of a Mat with a Mat
def  dot(that: Mat):  Mat = {
  var nv = new Mat(this.N, this.M )
   var i=0; var j=0;
   while (i<N) {
     j=0
    while (j<M)  {
      nv.v(i)(j) = v(i)(j) * that.v(i)(j)
      j += 1
    }
    i += 1
   }
 return nv
}

// dot product of a Mat with an Array[Array[Double]]
def  dot(that: Array[Array[Double]]):  Mat = {
  var nv = new Mat(this.N, this.M )
   var i=0; var j=0;
   while (i<N) {
     j=0
    while (j<M)  {
      nv.v(i)(j) = v(i)(j) * that(i)(j)
      j += 1
    }
    i += 1
   }
 return nv
}





// IN-PLACE Operations: Update directly the receiver, avoiding creating a new return object

    // Mat + Mat
def ++ (that: Mat): Mat =  {
  if (N != that.N || M != that.M)  // incompatible dimensions
      this
     else {
         var i=0; var j=0;
        while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) += that.v(i)(j)
      j +=1
    }
    i += 1
   }
   this
  }
}

    // Mat + Double
def ++ (that: Double): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) += that
      j +=1
    }
    i += 1
   }
   this
}

 // Mat + Int
def ++ (that: Int): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) += that
      j +=1
    }
    i += 1
   }
   this
}


    // Mat - Mat
def -- (that: Mat): Mat =  {
  if (N != that.N || M != that.M)  // incompatible dimensions
      this
     else {
         var i=0; var j=0;
        while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) -= that.v(i)(j)
      j +=1
    }
    i += 1
   }
   this
  }
}

    // Mat - Double
def -- (that: Double): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) -= that
      j +=1
    }
    i += 1
   }
   this
}

 // Mat - Int
def -- (that: Int): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) -= that
      j +=1
    }
    i += 1
   }
   this
}

     // Mat - Long
def -- (that: Long): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) -= that
      j +=1
    }
    i += 1
   }
   this
}


    // Mat * Double
def ** (that: Double): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) *= that
      j +=1
    }
    i += 1
   }
   this
}

    // Mat * Int
def ** (that: Int): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) *= that
      j +=1
    }
    i += 1
   }
   this
}


    // Mat * Long
def ** (that: Long): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) *= that
      j +=1
    }
    i += 1
   }
   this
}



    // Mat / Double
def /| (that: Double): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) /= that
      j +=1
    }
    i += 1
   }
   this
}

    // Mat / Int
def  /| (that: Int): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) /= that
      j +=1
    }
    i += 1
   }
   this
}


    // Mat / Long
def /| (that: Long): Mat =  {
  var i=0; var j=0
  while (i<N) {
       j=0
    while (j<M) {
      v(i)(j) /= that
      j +=1
    }
    i += 1
   }
   this
}




    // END OF IN-PLACE OPERATIONS


// Mat + Double
def + (that: Double): Mat =  {
      var nv = new Mat(this.N, this.M)
   var i=0; var j=0;
   while (i<N) {
       j=0
    while (j<M) {
      nv.v(i)(j) = v(i)(j)+that
      j +=1
    }
    i += 1
   }
   return nv
}

// Mat - Double
def - (that: Double): Mat =  {
      var nv = new Mat(this.N, this.M)
   var i=0; var j=0;
   while (i<N) {
     j=0
    while (j<M) {
      nv.v(i)(j) = v(i)(j)-that
      j += 1
    }
    i += 1
   }
   return nv
}


    // Mat + Double, for example at:   var m = ones0(3,4); var mm = 10.0 +: m,
  // the implicit conversion of 10.0 to RichNumber is avoided since the expression is evaluated as m.+(10.0)
def +: (that: Double): Mat =  {
      var nv = new Mat(this.N-1, this.M-1)
   var i=0; var j=0
   while (i<N) {
       j=0
    while (j<M) {
      nv.v(i)(j) = v(i)(j)+that
      j +=1
    }
    i += 1
   }
   return nv
}

    // Mat + Array[Double] : column oriented addition
def +: (that: Array[Double]): Mat =  {
    var vl = that.length
    if (vl != this.N) 
      return this

      var nv = new Mat(this.N-1, this.M-1)
   var i=0; var j=0
   while (j<M) {   // for all columns
       i=0
    while (i<N) {  // for all rows
      nv.v(i)(j) = v(i)(j)+that(i)
      j +=1
    }
    i += 1
   }
   return nv
}




    // Mat -< Vec : column oriented addition
def -< (that: Vec): Mat =  {
    var vl = that.length
    if (vl != this.N) 
      return this

      var nv = new Mat(this.N, this.M)
   var i=0; var j=0
   while (j<M) {   // for all columns
       i=0
    while (i<N) {  // for all rows 
      nv.v(i)(j) = v(i)(j)-that(i)
      i +=1
    }
    j += 1
   }
   return nv
}

   // Mat +< Vec : column oriented addition
def +< (that: Vec): Mat =  {
    var vl = that.length
    if (vl != this.N) 
      return this

      var nv = new Mat(this.N, this.M)
   var i=0; var j=0
   while (j<M) {   // for all columns
       i=0
    while (i<N) {  // for all rows 
      nv.v(i)(j) = v(i)(j)+that(i)
      i +=1
    }
    j += 1
   }
   return nv
}

   // Mat (NXM)*< Vec(MX1) : Mat(NX1) Matrix-Vector multiplication
def *< (that: Vec): Mat =  {
    var vl = that.length
    if (vl != this.M ) 
      return this
   

    var sm=0.0
    var nv = new Mat(this.N, 1)
   for (r<-0 to N-1) {      // all rows of the Mat
     sm=0.0
     for (c<-0 to M-1)
       sm += v(r)(c)*that(c)
     nv(r,0) = sm
    }
   
   nv
}

def -: (that: Double): Mat =  {
      var nv = new Mat(this.N-1, this.M-1)
   var i=0; var j=0
   while (i<N) {
       j=0
    while (j<M) {
      nv.v(i)(j) = v(i)(j)-that
      j +=1
    }
    i += 1
   }
   return nv
}


// Mat * Double
def *: (that: Double): Mat =  {
      var nv = new Mat(this.N-1, this.M-1)
   var i=0; var j=0
   while (i<N) {
     j=0
    while (j<M) {
      nv.v(i)(j) = v(i)(j)*that
      j += 1
    }
    i += 1
   }
   return nv
}


    // unary Minus applied to a Mat implies negation of all of its elements
def unary_- : Mat =  {
      var nv = new Mat(this.N, this.M)  // get a Mat of the same dimension
   var i=0; var j=0;
   while (i<N) {
     j=0
    while (j<M) {
      nv.v(i)(j) = -v(i)(j)  // negate element
      j += 1
    }
    i += 1
   }
   return nv
}

// transpose the Matrix
def trans: Mat = {
   var nv = new Mat(this.M, this.N)  // get a Matrix of dimension MXN
   var i=0; var j=0;
   while (i<N) {
     j=0
    while (j<M) {
      nv.v(j)(i) = v(i)(j)  // negate element
      j += 1
    }
    i += 1
   }

   return nv
}


// transpose the Matrix
def T: Mat = {
   var nv = new Mat(this.M, this.N)  // get a Matrix of dimension MXN
   var i=0; var j=0;
   while (i<N) {
     j=0
    while (j<M) {
      nv.v(j)(i) = v(i)(j)  // negate element
      j += 1
    }
    i += 1
   }

   return nv
}


    // multiply Mat*Double
def * (that: Double): Mat =  {
      var nv = new Mat(this.N, this.M)
   var i=0; var j=0;
   while (i<N) {
       j=0
    while (j<M)  {
      nv.v(i)(j) = v(i)(j)*that
      j+=1
    }
    i+=1
   }
   return nv
}

// Mat + Mat
def  + (that: Mat): Mat =  {
   var rN = this.N;   var rM = this.M;  // receiver's dimensions
   var sN = that.N;  var sM = that.M;  // parameter's dimensions

  if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Mat(sN, sM)
       var i=0; var j=0;
   while (i<sN)  {
        j=0
       while (j<sM)  {
           nv.v(i)(j) = v(i)(j) + that.v(i)(j)
           j+=1
     }
     i += 1
   }
   return nv
  }
  else { // incompatible dimensions
   that
     }
}




// Mat - Mat
def  -  (that: Mat): Mat =  {
   var rN = this.N;   var rM = this.M;
   var sN = that.N;  var sM = that.M;
   
    if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Mat(sN, sM)
       var i=0; var j=0;
   while (i<sN)  {
       j=0
    while (j<sM)  {
           nv.v(i)(j) = v(i)(j) - that.v(i)(j)
           j += 1
          }
         i += 1
    }
   return nv
  }  // same dimensions
  else { // incompatible dimensions
   return  that
     }
}

// Mat * Mat
 def * (that: Mat): Mat =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.N;  var sM = that.M;
   
   var  v1Colj = new Array[Double](rM)
   var result = new Mat(this.N, that.M)
   var j=0; var k=0;
   while (j < sM)  {
       k=0
      while  (k < rM) {
        v1Colj(k) =that.v(k)(j)
        k += 1
      }

      var i=0;
      while (i<rN) {
        var   Arowi = this.v(i)
        var   s = 0.0;
        k=0
        while (k< rM) {
          s += Arowi(k)*v1Colj(k)
          k += 1
        }
      result.v(i)(j) = s;
      i += 1
      }
 j += 1
   }
  return result
  }


// Mat + Array[Array[Double]]
 def + (that: Array[Array[Double]]): Mat =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.length;  var sM = that(0).length

    if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Mat(sN, sM)
       var i=0; var j=0;
   while (i<sN)  {
        j=0
       while (j<sM)  {
           nv.v(i)(j) = v(i)(j) + that(i)(j)
           j+=1
     }
     i += 1
   }
   return nv
  }
  else { // incompatible dimensions
   this
     }
}



// Mat - Array[Array[Double]]
 def - (that: Array[Array[Double]]): Mat =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.length;  var sM = that(0).length

 if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Mat(sN, sM)
       var i=0; var j=0;
   while (i<sN)  {
        j=0
       while (j<sM)  {
           nv.v(i)(j) = v(i)(j) - that(i)(j)
           j+=1
     }
     i += 1
   }
   return nv
  }
  else { // incompatible dimensions
   this
     }
}


// Mat * Array[Array[Double]]
 def * (that: Array[Array[Double]]): Mat =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.length;  var sM = that(0).length

  var  vr:Mat =new Mat(rN, sM)

   var   v1Colj = new Array[Double](rM)

   var j=0; var k=0;
   while (j < sM)  {
       k=0
      while  (k < rM) {
        v1Colj(k) =that(k)(j)
        k += 1
      }

      var i=0;
      while (i<rN) {
        var   Arowi = this.v(i)
        var   s = 0.0;
        k=0
        while (k< rM) {
          s += Arowi(k)*v1Colj(k)
          k += 1
        }
      vr.v(i)(j) = s;
      i += 1
      }
 j += 1
   }
  return vr
  }




//   the  M(fromR:incR:toR, fromC:incC:toC)  operation
def  subm(fromR:Int,  incR: Int,  toR: Int, fromC:Int,  incC: Int,  toC: Int)  =  
    apply(fromR, incR, toR, fromC, incC, toC)

//   the  M(fromR:incR:toR, :)  operation
def submr(fromR: Int, incR: Int, toR: Int) = 
    apply(':', fromR, incR, toR)

//   the  M(:, fromC:incC:toC)  operation
def submc(fromC: Int, incC: Int, toC: Int) = 
    apply(fromC, incC, toC, ':')



// postfix transpose operator
def ~() =  {
var  transposed = new Mat(this.M, this.N)
   for (r<-0 until this.M)
    for (c<-0 until  this.N)
      transposed(r,c) = this(c,r)
   transposed
} 


}


// Mat's companion object
  object Mat  {

   
 var mxElemsToDisplay = 6
 var matResizeFactor = 1.5
 var mxToStringElems = mxElemsToDisplay

// construct a Mat from a String , e,g,
// var m = M("3.4 -6.7; -1.2 5.6")
def  M(s: String) =  {
    var nRows = 1
    var nCols = 0
    for (i<-0 to s.length-1)   // count how many rows are specified
      if  (s(i)==';')
        nRows += 1

  // seperate rows to an ArrayBuffer
   var buf = new scala.collection.mutable.ArrayBuffer[String]()
   var strtok = new java.util.StringTokenizer(s, ";")  // rows are separated by ';'
   while (strtok.hasMoreTokens) {
         val tok = strtok.nextToken
         buf +=tok
      }    

// count how many numbers each row has. Assuming that each line has the same number of elements
 val firstLine = buf(0)
 strtok = new java.util.StringTokenizer(firstLine, ", ")  // elements are separated by ',' or ' '
 while (strtok.hasMoreTokens) {
   val tok = strtok.nextToken
   nCols += 1  
}

// allocate double array
var elems = new Array[Array[Double]](nRows, nCols)
  // pass elements from ArrayBuffer to double array
for (k <- 0 to buf.length-1)  {  // read array
   var currentLine = buf(k)
   strtok = new java.util.StringTokenizer(currentLine, ", ")  // elements are separated by ',' or ' '
var c=0 
while (strtok.hasMoreTokens) {  // read row
   val tok = strtok.nextToken
   var cNumber = tok.toDouble
   elems(k)(c) = cNumber
   c += 1
     }   // read row
   }  // read array

 // return the converted Matrix
 new Mat(elems)
 }  



// resamples the Mat every n rows and every m cols
def resample(matr: Mat, n: Int,  m: Int):Mat = {
    var matArr = matr.v
    var rows = matr.N
    var cols = matr.M
    var  rRows =  (rows/n).asInstanceOf[Int]
    var  rCols =  (cols/m).asInstanceOf[Int]
    var  newMatArr = new Array[Array[Double]](rRows, rCols)
    var r=0; var c=0
    while (r<rRows) {
        c = 0
        while  (c<rCols) {
            newMatArr(r)(c) = matr.v(r*n)(m*c)
            c += 1
          }
      r += 1
}
   var   resampledMat = new Mat(newMatArr);
    resampledMat
}

// convert Vec to Mat
def vecToMat( that: Vec): Mat = {   
  println("converting Vec")
  var N = that.length
  var rmat = new Mat(N, 1)
  for (k<-0 to N-1)
   rmat(k, 0) = that(k)

 rmat
}

  // COLUMNWISE OPERATIONS
// columnwise sum
def sum(matr: Mat): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     sm=0.0
     for (crow<- 0 to N-1)
       sm += matr.v(crow)(ccol)
     res(ccol) = sm
     }
 res
}

// columnwise mean
def mean(matr: Mat): Array[Double]  = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](M)
    for (ccol <-0  to M-1) {
     sm=0.0
     for (crow <- 0 to N-1)
       sm += matr.v(crow)(ccol)
     res(ccol) = sm/N
     }
 res
}

// columnwise product
def prod(matr: Mat): Array[Double]  = {
    var N = matr.N;     var M = matr.M
    var pd = 1.0
    var res = new Array[Double](M)
    for (ccol <- 0  to M-1) {
     pd=1.0
     for (crow <- 0 to N-1)
       pd *= matr.v(crow)(ccol)
     res(ccol) = pd
     }
    res
}

 // columnwise min
 def min(matr: Mat): Array[Double] = {
     var N = matr.N;   var M = matr.M
     var res = new Array[Double](M)
     var ccol = 0
     while  (ccol <= M-1) {
         var mn = matr.v(0)(ccol)  // keeps the running min element
         var crow = 0
         while  (crow <= N-1)  {
              var tmp = matr.v(crow)(ccol)
               if (tmp  < mn)  mn = tmp
               crow += 1
           }
           res(ccol) = mn   // min element for the ccol column
           ccol += 1
     }
     res
 }

  // columnwise max
 def max(matr: Mat): Array[Double] = {
     var N = matr.N;   var M = matr.M
     var res = new Array[Double](M)
     for (ccol <- 0 to M-1) {
         var mx = matr.v(0)(ccol)  // keeps the running max element
         for (crow <- 1 to N-1)  {
              var tmp = matr.v(crow)(ccol)
               if (tmp  > mx)  mx = tmp
           }
        res(ccol) = mx   // max element for the ccol column
     }
     res
 }


// ROWWISE OPERATIONS
// rowwise sum
def sumR(matr: Mat): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](N)  // sum for all rows
    for (crow<- 0 to N-1)  {
        sm=0.0
        for (ccol <-0  to M-1)   // sum across column
         sm += matr.v(crow)(ccol)
     res(crow) = sm
     }
 res
}

// rowwise mean
def meanR(matr: Mat): Array[Double]  = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](N)  // mean for all rows
    for (crow <-0  to N-1) {
     sm=0.0
     for (ccol<- 0 to M-1)  // sum across column
       sm += matr.v(crow)(ccol)
     res(crow) = sm/M
     }
 res
}

// rowwise product
def prodR(matr: Mat): Array[Double]  = {
    var N = matr.N;     var M = matr.M
    var pd = 1.0
    var res = new Array[Double](N)
    for (crow <-0  to N-1) {
     pd=1.0
     for (ccol<- 0 to M-1)  // product across column
       pd *= matr.v(crow)(ccol)

     res(crow) = pd
     }
    res
}

 // rowwise min
 def minR(matr: Mat): Array[Double] = {
     var N = matr.N;   var M = matr.M
     var res = new Array[Double](N)
     for (crow <- 0 to N-1) {
         var mn = matr.v(crow)(0)  // keeps the running min element
         for (ccol <- 1 to M-1)  {
              var tmp = matr.v(crow)(ccol)
               if (tmp  < mn)  mn = tmp
           }
           res(crow) = mn   // min element for the crow row
     }
      res
 }

  // rowwise max
 def maxR(matr: Mat): Array[Double] = {
     var N = matr.N;   var M = matr.M
     var res = new Array[Double](N)
     for (crow <- 0 to N-1) {
         var mx = matr.v(crow)(0)  // keeps the running max element
         for (ccol <- 1 to M-1)  {
              var tmp = matr.v(crow)(ccol)
               if (tmp  > mx)  mx = tmp
           }
        res(crow) = mx   // max element for the ccol column
     }
     res
 }



   def diag( x: Array[Double]) = {
    var n = x.length
    var om = new Mat(n, n)
    var i=0;
    while  (i< n) {
      om.v(i)(i) = x(n)
       i += 1
      }
          om
    }

     def diag(n: Int): Mat = {
         var om = new Mat(n, n)
         var i=0;
         while  (i< n) {
               om.v(i)(i) = 1.0
               i += 1
              }
          om
    }

     def ones(n: Int, m:Int): Mat = {
         var om = new Mat(n, m)
         var i=0; var j=0;
         while  (i< n) {
            j=0
            while (j < m) {
               om.v(i)(j) = 1.0
               j += 1
              }
            i += 1
         }
          om
    }

      def zeros(n: Int, m:Int): Mat = {
         var om = new Mat(n, m)
         var i=0; var j=0;
         while  (i< n) {
            j=0
            while (j < m) {
               om.v(i)(j) = 0.0
               j += 1
              }
            i += 1
         }
          om
    }

 def fill(n: Int, m:Int, vl: Double): Mat = {
         var om = new Mat(n, m)
         var i=0; var j=0;
         while  (i< n) {
            j=0
            while (j < m) {
               om.v(i)(j) = vl
               j += 1
              }
            i += 1
         }
          om
    }

 
def  eye(n:Int): Mat = {
    eye(n, n)
}

def  eye(n:Int, m:Int): Mat = {
     var om = new Mat(n, m)
     var minCoord = n
     if (m<minCoord)
        minCoord = m
     var i=1 
     while  (i< minCoord) {
             om.v(i)(i) = 1.0
              i += 1
         }
          om
     }



    def rand(n: Int, m:Int): Mat = {
        var om = new Mat(n, m)
       var i=0; var j=0;
       while  (i< n) {
            j=0
            while (j < m) {
             om.v(i)(j) = java.lang.Math.random
             j += 1
           }
          i += 1
       }
          om

    }



 def sin(v: Mat): Mat = {
       var N = v.N; var M = v.M;
        var om = new Mat(N, M)
       var i=0; var j=0;
       while  (i< N) {
            j=0
            while (j < M) {
               om.v(i)(j) = java.lang.Math.sin(v.v(i)(j))
               j += 1
            }
            i += 1
       }
          om
    }

 def sin(vin: Array[Double]): Array[Double] = {
      var  n = vin.length
      var  v = new Array[Double](n)
      var i=0
      while (i<n)  {
          v(i) =  java.lang.Math.sin(vin(i))
          i += 1
      }
      v
  }



    def cos(v: Mat): Mat = {
        var N = v.N; var M = v.M;
        var om = new Mat(N, M)
        var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
             om.v(i)(j) = java.lang.Math.cos(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def tan(v: Mat): Mat = {
      var N = v.N; var M = v.M;
        var om = new Mat(N, M)
        var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
             om.v(i)(j) = java.lang.Math.tan(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



def asin(v: Mat): Mat = {
      var N = v.N; var M = v.M;
        var om = new Mat(N, M)
        var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
             om.v(i)(j) = java.lang.Math.asin(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def asin(v: Array[Array[Double]]): Array[Array[Double]] = {
  var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.asin(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


def acos(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
           j=0
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.acos(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

 def acos(v: Array[Array[Double]]): Array[Array[Double]] = {
   var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.acos(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



def atan(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
           j=0
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.atan(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def   atan(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.atan(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


def sinh(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
           j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.sinh(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


def sinh(v: Array[Array[Double]]): Array[Array[Double]] = {
      var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.sinh(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }

def cosh(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
          j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.cosh(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



def cosh(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.cosh(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }


def tanh(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
          j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.tanh(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
 }

def tanh(v: Array[Array[Double]]): Array[Array[Double]] = {
       var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
             j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.tanh(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }



def exp(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
           j=0
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.exp(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def exp(v: Array[Array[Double]]): Array[Array[Double]] = {
      var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.exp(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }

def log(v: Mat): Mat = {
  var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
          j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.log(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

   def log(v: Array[Array[Double]]): Array[Array[Double]] = {
        var N = v.length; var M = v(0).length;
        var  vd =new Array[Array[Double]](N)
        var c = 0
        while (c < N) {
         vd(c) = new Array[Double](M)
         c += 1
        }
       var i=0; var j=0;
        while (i < N) {
            j=0
            while (j < M)  {
              vd(i)(j) = java.lang.Math.log(v(i)(j))
              j += 1
            }
            i += 1
        }
          vd
    }

 def log2(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
var conv = java.lang.Math.log(2.0)
   var i=0; var j=0;
     while (i < N) {
           j=0
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.log(v.v(i)(j)/conv)
             j += 1
            }
            i += 1
        }
          om
    }


  

def  log10(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
         j=0
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.log10(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }




def abs(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
          j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.abs(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


 def ceil(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
         j=0
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.ceil(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def floor(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
          j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.floor(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



 def round(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
         j=0
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.round(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


def sqrt(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
          j=0
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.sqrt(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



 def toDegrees(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
           j=0
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.toDegrees(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



 def toRadians(v: Mat): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
         j=0
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.toRadians(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



 def pow(v: Mat, exponent:Double): Mat = {
   var N = v.N; var M = v.M;
   var om = new Mat(N, M)
   var i=0; var j=0;
     while (i < N) {
           j=0
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.pow(v.v(i)(j), exponent)
             j += 1
            }
            i += 1
        }
          om
    }





 def  randomNormal(m: Int,  n: Int,  mu:Double, sigma:Double): Mat  = {
		var g =  scalaSci.math.array.StatisticSample.randomNormal(m, n, mu, sigma)
        var matRN = new Mat(g)
        matRN
}

def  randomUniform(m: Int,  n: Int,  min:Double, max:Double): Mat  = {
		var g =  scalaSci.math.array.StatisticSample.randomUniform(m, n, min, max)
        var matRN = new Mat(g)
        matRN
}

////////   scalaSci.math.LinearAlgebra.LinearAlgebra routines
def CholeskyL(Mt: Mat): Mat = {
 var choleskDec = LinearAlgebra.cholesky(Mt.v)
 var choleskMatDoubles = choleskDec.getL().getArray()
 new Mat(choleskMatDoubles)
}

def  Cholesky_SPD(Md: Mat): Mat = {
      var CholeskyDecomposition =  LinearAlgebra.cholesky(Md.v)
      var  choleskMatDoubles = CholeskyDecomposition.getL().getArray()
      new Mat(choleskMatDoubles)
}

def  Cholesky_solve(Md: Mat,  b: Mat ) = {
      var   choleskDec = LinearAlgebra.cholesky(Md.v)
      var   jb = new jMatrix(b.v)
      var   solvedMat = choleskDec.solve(jb)
      new Mat(solvedMat.getArray())
}

def  LU_L(M: Mat):Mat = {
    new Mat(LinearAlgebra.LU( M.v).getL().getArray())
}

def  L(M: Mat): Mat = {
   LU_L( M)
}


def  LU_U(M: Mat): Mat = {
   new Mat(LinearAlgebra.LU( M.v).getU().getArray())
}

def U(M: Mat): Mat = {
   LU_U( M)
 }

def LU_det(M: Mat) = {
   LinearAlgebra.LU( M.v).det()
}

def LU(A: Mat):(Mat, Mat, Mat) = { 
  val LUdecomp = LinearAlgebra.LU( A.v)   // perform the LU decomposition
  var L =  new Mat(LUdecomp.getL().getArray())
  var U =  new Mat(LUdecomp.getU().getArray())
  var P = new Mat(A.N, A.N)
  var p = LUdecomp.getPivot()
 // construct permutation matrix
    for (k<-0 until A.N)  
      P(k, p(k)) = 1
 
  (L, U, P)
 } 


def  det(M: Mat) = {
  LinearAlgebra.det(M.v)
}

def  LU_solve(A:Mat,  b: Mat): Mat = {
    var  LUDec = LinearAlgebra.LU(A.v)
    var  jb = new jMatrix(b.v)
    var  solvedMat = LUDec.solve(jb)
    new Mat(solvedMat.getArray())
}

def  QR_H(M: Mat): Mat = {
   new Mat(LinearAlgebra.QR( M.v).getH())
}

def  QR_Q(M: Mat): Mat = {
  new Mat(LinearAlgebra.QR( M.v).getQ())
}

def  Q(M: Mat): Mat = {
    QR_Q( M)
 }

def   QR_R(M: Mat) = {
  new Mat(LinearAlgebra.QR( M.v).getR())
}

def  R(M: Mat) =  {
   QR_R( M)
 }

def QR(M: Mat) = {
  (Q(M), R(M))
}

def  QR_solve(A: Mat, b: Mat) = {
  var  QRDec = LinearAlgebra.QR(A.v)
  var  jb = new jMatrix(b.v)
  var  solvedMat = QRDec.solve(jb.getArray())
  new Mat(solvedMat)
}

def solve(A: Mat,  b: Mat) = {
   new Mat(LinearAlgebra.solve( A.v, b.v))
}

def Singular_cond(M: Mat): Double = {
   LinearAlgebra.singular( M.v).cond()
}

def cond(M: Mat): Double = {
   LinearAlgebra.cond( M.v)
}

def  Singular_S(M: Mat): Mat = {
   new Mat(LinearAlgebra.singular(M.v).getS().getArray())
}

def  S(M: Mat): Mat = {
  Singular_S( M)
}

def  Singular_values(M: Mat): Mat = {
  new Mat(LinearAlgebra.singular(M.v).getSingularValues())
}

def  Singular_U(M: Mat): Mat = {
  new Mat(LinearAlgebra.singular(M.v).getU().getArray())
}

def   Singular_V(M: Mat): Mat = {
  new Mat(LinearAlgebra.singular(M.v).getV().getArray())
}

def Singular_norm2(M: Mat): Double = {
  LinearAlgebra.singular( M.v).norm2()
}

def norm2(M: Mat): Double = {
  LinearAlgebra.norm2( M.v)
}

def  norm1(M: Mat): Double = {
  LinearAlgebra.norm1( M.v)
}

def  normF(M: Mat): Double = {
  LinearAlgebra.normF( M.v)
}

def  normInf(M: Mat): Double = {
   LinearAlgebra.normInf( M.v)
}

def  Singular_rank(M: Mat): Int = {
   LinearAlgebra.singular( M.v).rank()
}

//    var ([U,S,V) = svd(X) produces a diagonal matrix S, of the same 
//    dimension as X and with nonnegative diagonal elements in
//     decreasing order, and unitary matrices U and V so that    X = U*S*V'.
 def svd(X: Mat) = {
  // perform the Singular Value Decomposition 
   
  val svdM = new Jama.SingularValueDecomposition(new Jama.jMatrix(X.getv))
  val U = svdM.getU.getArray
  val V = svdM.getV.getArray
  val s = svdM.getSingularValues
  val N = s.length
  val S = new Array[Array[Double]](N, N)
  for (k<-0 until N)
     S(k)(k) = s(k)
  (U, S, V)

}


def  rank(M: Mat): Int = {
  LinearAlgebra.rank( M.v)
}

def trace(M: Mat): Double = {
   LinearAlgebra.trace( M.v)
}

def Eigen_V(M: Mat): Mat = {
   new Mat(LinearAlgebra.eigen(M.v).getV())
}

def V(M: Mat): Mat = {
   Eigen_V( M)
 }

def  Eigen_D(M: Mat) = {
   new Mat(LinearAlgebra.eigen(M.v).getD())
}

def  D(M: Mat):Mat = {
   Eigen_D( M)
}
// return eigenvectors, eigenvalues
def eig(M: Mat) = {
 (V(M), D(M))
}

def inv(M: Mat): Mat = {
    var jm = new Jama.jMatrix (M.v)
    var invm = jm.inverse
    new Mat(invm.getArray)
}

 def cov(mt: Mat): Mat = {
     var v = mt.v
     var r = scalaSci.math.array.StatisticSample.covariance(v, v);
     new Mat(r)
 }

 def covariance(mt: Mat): Mat = {
     var v = mt.v
     var r = scalaSci.math.array.StatisticSample.covariance(v, v);
     new Mat(r)
 }



 def  corr(mt: Mat): Mat = {
     var v =  mt.v
     var r = scalaSci.math.array.StatisticSample.correlation(v, v);
     new Mat(r)
 }


def  transpose(Mt: Mat): Mat =  {
    return new Mat(DoubleArray.transpose(Mt.v ) )
	}

  def  T(Mt: Mat): Mat =  {
    return new Mat(DoubleArray.transpose( Mt.v ) )
	}


// construct a Mat from a double [][] array, one indexed
def MatFromOneIndexedDArr(vals: Array[Array[Double]]): Mat = {
       var vals0 = new Array[Array[Double]](vals.length-1, vals(0).length-1)
    var i=1; var j=1;
    while (i<vals.length) {
        j=1
        while (j<vals(0).length) {
       vals0(i-1)(j-1) = vals(i)(j)
       j+=1
        }
        i+=1
     }
     new Mat(vals0)
 }



//  Reshapes a  matrix a to new dimension n X m
def  reshape(a: Mat,  n: Int, m: Int): Mat = {
   var  aCols  = a.M   // columns of matrix
   var  aRows = a.N   // rows of matrix
   var N = n; var M = m;
  // m,n must be positive
 if (M<=0 || N<=0)  {  N = 1; M = aRows*aCols;  }  // default to reshaping in a large row
 if (N*M != aRows*aCols)  {
     return   a  // invalid new size: return the original matrix
    }
 var  nm = new Mat(N, M)  // create the new matrix
// keep two set of indices, i.e. iorig, jorig: indices at the original matrix and inew, jnew: indices at the new matrix
var iorig=0; var jorig=0; var inew=0; var jnew=0
while (iorig < aRows)  {
    while (jorig < aCols)  {
        nm(inew, jnew) = a(iorig, jorig)
        jorig += 1
        jnew += 1
        if (jnew>=M) {  // next row of the reshaped new matrix
            jnew = 0
            inew += 1
           }
      } // jorig < aCols
    iorig += 1
    jorig = 0
  } // iorig < aRows
nm  // return the new Matrix
}

def testMat(N: Int, M: Int) = {
  var a = new Mat(N, M)
  for (rows <- 0 until N)
    for (cols <-0 until M)
      a(rows, cols) = rows*10+cols

  a
}

  

}

